Start date: 04-15-2020
End date: 04-16-2020
Analysed by: Ruth Chia
Working directory: /data/NDRS_LNG/MyastheniaGravis/updated.April2020/Itals
Run sample- and variant-level QC on raw genotype plink files.
Submit for imputation using Topmed Imputation Server (hg38).
Generate PCs and covariate file.
Run GLM (dosages) using imputed data
QC script is written in two parts i.e.
Sample-level QC includes:
STEP 1 : SAMPLE-LEVEL FILTERING (HETEROZYGOSITY, CALL RATE, GENDER CHECK, ANCESTRY)
STEP 1.1: Sample-level filtering (Heterozygosity check, removal of het outliers)
STEP 1.2: Sample-level filtering (sample call rate check, maximum missing calls per sample --mind 0.05)
STEP 1.3: Sample-level filtering (Gender check, remove samples that failed gender check)
STEP 1.4: Sample-level filtering (Ancestry check, remove samples that are not european)
STEP 1.5: Generate list of duplicates; remove duplicates
STEP 1.6: Generate list of related samples
Variant-level QC includes:
STEP 2 : VARIANT-LEVEL FILTERING (Variant inclusion criteria to account for genotyping batch differences; then make final bfiles with --geno 0.05)
STEP 2.1: Case/control nonrandom missingness test (i.e. exclude SNPs with P \< 1E-4)
STEP 2.2: Haplotype-based test for non-random missing genotype data (i.e. exclude SNPs with P \< 1E-4)
STEP 2.3: Make list of varaints that failed HWE (controls) (i.e. exclude SNPs with midP \< 1E-10)
STEP 2.4: Final --geno 0.05 to account for final variant missingness
%%bash
echo "sh ./scripts/QC_preimputation_v2_keepRelated_part1of2.mg.sh merged1 Itals_mg chiarp@mail.nih.gov" > qc1.swarm
swarm --file qc1.swarm --gres=lscratch:200 --logdir swarmOE_qc --partition quick \
-g 64 -t auto --time 4:00:00 --sbatch "--mail-type=BEGIN,FAIL,TIME_LIMIT_80"
move QC-ed files to specific directory
%%bash
module load plink/1.9.0-beta4.4
mkdir CLEAN.rawGenotype.keepRelated
mv FILTERED.Itals_mg_keepRelated.* ./CLEAN.rawGenotype.keepRelated
cd CLEAN.rawGenotype.keepRelated
mv FILTERED.Itals_mg_keepRelated.bed FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005.bed
mv FILTERED.Itals_mg_keepRelated.bim FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005.bim
mv FILTERED.Itals_mg_keepRelated.fam FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005.fam
%%bash
module load plink/1.9.0-beta4.4
mkdir CLEAN.rawGenotype.UNRELATED
plink \
--bfile CLEAN.rawGenotype.keepRelated/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005 \
--remove FILTERED.to.remove.GRM_matrix.noDups.relatedSamples_FIDspaceIID.txt \
--keep-allele-order \
--allow-no-sex \
--make-bed \
--out CLEAN.rawGenotype.UNRELATED/FILTERED.Itals_mg_noDups.UNRELATED.hwe1e-10.geno005
%%bash
cp CLEAN.rawGenotype.keepRelated/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005.* .
echo "./scripts/QC_preimputation_v2_keepRelated_part2of2.sh FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005 chiarp@mail.nih.gov" > QC_preimputation_part2.swarm
swarm --file QC_preimputation_part2.swarm --gres=lscratch:200 \
--module plink,R -g 32 --time 24:00:00 -t auto \
--sbatch "--mail-type=BEGIN,FAIL,END,TIME_LIMIT_80"
%%bash
# Remove some intermediate files that are not needed to save some space
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr*.check.*
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005X.*
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005[0-9].*
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005[0-9][0-9].*
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-updated-chr*.*
rm *-FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-HRC.txt
rm Run-plink.sh
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005.*
To create covariate files for:
RELATED+UNRELATED samples UNRELATED samples only%%bash
module load plink/1.9.0-beta4.4
module load R/3.5.2
module load GCTA/1.92.0beta3
module load flashpca/2.0
mkdir PCAcovs
# Prune snps
plink \
--bfile CLEAN.rawGenotype.keepRelated/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005 \
--allow-no-sex \
--maf 0.01 \
--geno 0.01 \
--hwe 5e-6 \
--autosome \
--exclude range ./scripts/exclusion_regions_hg19.txt \
--make-bed \
--out PCAcovs/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter \
--memory 119500 --threads 19
plink \
--bfile PCAcovs/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter \
--allow-no-sex \
--geno 0.01 \
--maf 0.05 \
--indep-pairwise 1000 10 0.02 \
--out PCAcovs/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter_pruning \
--memory 119500 --threads 19
# All samples (related + unrelated samples; no dups)
plink \
--bfile PCAcovs/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter \
--allow-no-sex \
--extract PCAcovs/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter_pruning.prune.in \
--keep-allele-order \
--make-bed \
--out PCAcovs/pruned.Itals_mg_noDups.keepRelated \
--memory 119500 --threads 19
# Unrelated samples only (no dups)
plink \
--bfile PCAcovs/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter \
--remove GRM/FILTERED.to.remove.GRM_matrix.noDups.relatedSamples_FIDspaceIID.txt \
--allow-no-sex \
--extract PCAcovs/FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter_pruning.prune.in \
--keep-allele-order \
--make-bed \
--out PCAcovs/pruned.Itals_mg_noDups.UNRELATED \
--memory 119500 --threads 19;
# Calculate/generate PCs based on pruned data set
cd PCAcovs
flashpca --bfile pruned.Itals_mg_noDups.keepRelated -d 10 --suffix _Itals_mg_noDups.keepRelated_forPCA --numthreads 19
flashpca --bfile pruned.Itals_mg_noDups.UNRELATED -d 10 --suffix _Itals_mg_noDups.UNRELATED_forPCA --numthreads 19
# Move all log files to a new folder
mkdir logFiles
mv *.log logFiles
# Remove intermediate files
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter.bed
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter.bim
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter.fam
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter_pruning.prune.in
rm FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005_filter_pruning.prune.out
Create covariate file
Use updated phenotype file as input for here: phenotype_mg_Carlo2018_2012_ItalsControls.updatedMarch2020.FINAL.txt. Then subset to only keep samples that made it to the 'CLEAN' plink1 files.
%%bash
module load R/3.5.2
R --vanilla --no-save
# Load libraries
require(data.table)
require(tidyverse)
# UNRELATED samples only (no dups)
## Read in files
pc <- fread("PCAcovs/pcs_Itals_mg_noDups.UNRELATED_forPCA", header=T)
fam <- fread("PCAcovs/pruned.Itals_mg_noDups.UNRELATED.fam", header=F)
fam <- fam %>% rename(FID = V1, IID = V2, MAT = V3, PAT = V4, SEX = V5, PHENO= V6)
pheno <- fread("phenotype_mg_Carlo2018_2012_ItalsControls.updatedMarch2020.FINAL.txt", header=T)
pheno <- pheno %>% rename(genderONFILE = gender)
dim(pc)
dim(fam)
dim(pheno)
## Merge pc, fam and pheno files.
## Also create additional column so that the sample ID matches the format in the dose.vcf.gz (i.e. FID_IID)
data1 <- merge(fam,pheno,by=c("FID","IID"))
data2 <- merge(data1, pc, by=c("FID","IID"))
data2$GENDER <- data2$SEX - 1
data2$genderONFILE[is.na(data2$genderONFILE)] <- "NA"
data2$FID_IID <- paste(data2$FID,data2$IID, sep="_")
data2$age_at_onset <- as.numeric(data2$age_at_onset)
write.table(data2, "COVARIATES.Itals_mg_noDups.UNRELATED.txt", sep="\t",quote=F, col.names=T,row.names=F)
data3 <- data2 %>%
mutate(FID = FID_IID, IID = FID_IID) %>%
select(-FID_IID)
write.table(data3, "COVARIATES.Itals_mg_noDups.UNRELATED.forGLM.txt", sep="\t",quote=F, col.names=T,row.names=F)
head(data3)
## Summarise numbers by pheno or diagnosis and gender or sex
## Note that there is a mismatch of gender from pheno file compared to the actual SEX information in the plink file.
## But since these samples survived gender check.
## Decided to stick with information in the .fam file
table(data2$PHENO, data2$SEX)
table(data2$diagnosis, data2$SEX)
table(data2$diagnosis, data2$genderONFILE)
table(data2$diagnosis, data2$GENDER)
head(data2)
## Run step() to determine which covariate to correct for for association analysis
data2$CASE <- data2$PHENO - 1
unrelated <- glm(CASE ~ GENDER + age_at_onset + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 , data = data2, family = "binomial"(link = "logit"))
summary(unrelated)
step(unrelated)
UNRELATED samples only (no dups):
3567.27GENDER,age_at_onset,PC1,PC2,PC3,PC4,PC5,PC6,PC10%%bash
# Plot PCs to check to make sure cases and controls are not clustering out
module load R/3.5.2
Rscript ./scripts/PCplots_MG.GWAS.R COVARIATES.Itals_mg_noDups.keepRelated.txt PCplots_Itals_mg_noDups.keepRelated
Rscript ./scripts/PCplots_MG.GWAS.R COVARIATES.Itals_mg_noDups.UNRELATED.txt PCplots_Itals_mg_noDups.UNRELATED
from IPython.display import display
from PIL import Image
print("PC plots for samples and cases (UNRELATED samples only)")
PCplot="PCplots_Itals_mg_noDups.UNRELATED_prunedSNPs_hwe5e-6.geno001.maf001.jpg"
display(Image.open(PCplot))
Uploaded .vcf.gz files for imputation:
pre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr1.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr2.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr3.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr4.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr5.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr6.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr7.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr8.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr9.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr10.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr11.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr12.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr13.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr14.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr15.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr16.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr17.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr18.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr19.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr20.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr21.vcf.gzpre_impute_FILTERED.Itals_mg_noDups.keepRelated.hwe1e-10.geno005-chr22.vcf.gz%%bash
echo "sh scripts/MG.ITALS_ImputationResultsDownload.hg38_and_plink_04-28-2020.sh" > MG.ITALS_ImputationDownload.hg38.swarm
swarm --file MG.ITALS_ImputationDownload.hg38.swarm \
--logdir swarmOE_ImputationDownload \
--gres=lscratch:200 \
--module plink/1.9.0-beta4.4 \
-g 2 -t auto --time 24:00:00 \
--sbatch "--mail-type=BEGIN,FAIL,TIME_LIMIT_80"
With the output from minimac 4, unzipped files does not automatically create an index file. Do this manually.
%%bash
cd Imputation.hg38
for CHNUM in {1..22};
do
echo "tabix -p vcf chr${CHNUM}.dose.vcf.gz" >> tabix_dose.vcf.gz.swarm
done
swarm --file tabix_dose.vcf.gz.swarm --gres=lscratch:400 \
--logdir swarmOE_tabix_dose.vcf.gz \
--module samtools -g 120 -t auto --time 72:00:00 \
--sbatch "--mail-type=BEGIN,FAIL,TIME_LIMIT_80"
%%bash
DIR="/data/NDRS_LNG/MyastheniaGravis/updated.April2020/Itals"
cd $DIR/Imputation.hg38
wc -l chr[0-9]*.info
%%bash
mkdir Imputation.hg38/Vars.Rsq03MAF00001
cd Imputation.hg38
for CHNUM in {1..22}
do
awk '{if($5 > 0.0001 && $7 > 0.3) print $1}' chr${CHNUM}.info | tail -n +2 > Vars.Rsq03MAF00001/Rsq03MAF00001.chr${CHNUM}.txt
done
!wc -l Imputation.hg38/Vars.Rsq03MAF00001/Rsq03MAF00001.chr[0-9]*.txt
What is needed:
Location of files:
USmerged
/data/NDRS_LNG/MyastheniaGravis/updated.April2020/US/Analysis.GLM.hg38/US.JointPostImputation/merged.vcf/US.chr${CHNUM}.vcf.gz/data/NDRS_LNG/MyastheniaGravis/updated.April2020/US/Analysis.GLM.hg38/US.JointPostImputation/sharedVars.postImputation/SharedVars.postImputation.Rsq03MAF00001.chr${CHNUM}.txt/data/NDRS_LNG/MyastheniaGravis/updated.April2020/US/COVARIATES.US_mg_noDups.UNRELATED.forImputed.txtItaly
/data/NDRS_LNG/MyastheniaGravis/updated.April2020/Itals/Imputation.hg38 /chr${CHNUM}.dose.vcf.gz/data/NDRS_LNG/MyastheniaGravis/updated.April2020/Itals/Imputation.hg38/Vars.Rsq03MAF00001/Rsq03MAF00001.chr${CHNUM}.txt/data/NDRS_LNG/MyastheniaGravis/updated.April2020/Itals/COVARIATES.Itals_mg_noDups.UNRELATED.forGLM.txt%%bash
cd Imputation.hg38
mkdir plinkForPRS
DIR="/data/NDRS_LNG/MyastheniaGravis/updated.April2020/Itals"
# Make list of unrelated samples to keep
awk '{print $1,$2}' $DIR/COVARIATES.Itals_mg_noDups.UNRELATED.forGLM.txt > $DIR/Imputation.hg38/SampleList.Itals.UNRELATED.forImputed.txt
for CHNUM in {1..22}
do
echo "sh vcfImputedToplink1.Itals.sh $DIR $CHNUM" >> generatePRSplink.swarm
done
swarm --file generatePRSplink.swarm --logdir swarmOE_vcfToPlink -g 120 \
--time 2:00:00 -t auto --module plink/1.9.0-beta4.4 --partition quick
%%bash
# view log file
cd Imputation.hg38/
cat swarmOE_vcfToPlink/swarm_10041423_0.o